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ABSTRACT 



In Augvist 1992, an acoustic tomography experiment was conducted in the 
Barents Sea. The objective of this thesis is the development of space-time signal processing 
algorithms to extract the arrival times of rays and modes for that experiment. The 
temporal signal processing involves pulse compression of the phase-encoded signal. 
Spatial signal processing consists of plane wave and modal beamforming on data from a 
16 element vertical line array. After signal processing an arrival detection, identification, 
and tracking algorithm is used to obtain travel time and travel time error estimates. After 
development and testing of these algorithms, 13 hours of acoustic data was processed for 
rays and eight hours of data was processed for modes. The resultant travel time 
information is then available for inversion to sound speed maps. 
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I. INTRODUCTION 



A. ACOUSTIC TOMOGRAPHY IN THE LITTORAL REGION 

Future conflicts will require naval forces to operate in coastal regions [Ref. 1]. In 
order to achieve battle space dominance in these regions, the subsurface environment must 
be secured first. Securing the subsurface environment will require success in the 
prosecution of diesel submarines and in the hunting of mines. A key to success will be 
a detailed knowledge of ocean acoustic propagation. Traditional methods of obtaining the 
necessary environmental information will not provide the spatial resolution needed in a 
tactical time frame. Traditional methods can not provide sonar performance prediction 
of sufficient quality. Acoustic tomography provides a method to obtain the environmental 
information needed to determine and/or predict acoustic propagation with the necessary 
spatial resolution in a tactical time frame for the shallow water environment. Timely high 
resolution sound speed maps will be needed to achieve battle space dominance in littoral 
regions. 

B. MOVING TOMOGRAPHY FROM DEEP TO SHALLOW WATER 

Ocean acoustic tomography was originally proposed by Munk and Wunsch in 1977. 

In 1979, they presented methods for inverting the data to estimate the sound speed field 
[Ref. 2]. The procedure is similar to that used in medical x-ray tomography where the 
measuring signal travels in a straight line from source to receiver. Ocean acoustic 
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tomography relies on acoustic energy traveling over many paths with different travel 
times. 

The word tomography is derived from two Greek roots meaning "to slice" and "to 
look at". Acoustic tomography is then looking at a slice of the ocean with acoustics. 
Sound speed is a function of salinity, pressure, and temperature. As acoustic energy 
travels along different paths, its speed depends upon these physical quantities. Travel time 
along these paths can be used to find a spatial sound speed map of the ocean. The 
knowledge of how sound speed affects travel time is known as i\\t forward problem. The 
process of estimating sound speed from the travel time is known as the inverse problem 
or inversion [Ref. 2]. 

Tomography has been successfully used in deep water using receivers consisting of 
a single hydrophone or very simple arrays. Deep water tomography has been used to map 
regions of the deep ocean over ranges of 300 to 1000 km. This long travel distance 
provides sufficient separation in time for multipath arrival identification. Ray models are 
able to accurately describe much of the propagation. In shallow water, the acoustic 
propagation becomes more complex making the forward and inverse problems more 
difficult. Surface and bottom interactions occur frequently and multipath arrival 
identification and time estimation is made difficult by the lack of temporal separation. 
To identify these multipath arrivals, a vertical array can be used to provide temporal and 
spatial information on these arrivals. A vertical array can also be used as a sparial filter 
to beamform for the acoustic normal modes [Ref 3, 4]. Thus a modal model can be used 
for propagation and inversion in addition to the ray model. 
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The primary advantage of acoustic tomography is its ability to sample an ocean 
volume quickly, i.e., at the speed of sound. For a typical shallow water survey, it may 
take many hours with traditional oceanographic sampling techniques to determine a sound 
speed field. Once the acoustic tomography source and receiver are in place, sound energy 
takes seconds to travel from the source to the vertical array. A sound speed field can be 
calculated as often as desired subject to the battery limitations of the sources and 
receivers. Tomography provides a unique method to view the interior of the ocean and 
its variability. Satellite remote sensing systems may be complimentary to tomography by 
providing information on features visible at the surface where resolution of the 
tomography maps may be poor. In addition, both tomography and satellite remote 
sensing data can be assimilated into an ocean model potentially providing very high 
resolution maps of a coastal ocean volume [Ref 3]. 

C. OVERVIEW OF THESIS 

The tomography problem is summarized in Figure 1.1. As Figure 1.1 shows, 
imaging the ocean sound speed variability with acoustic tomography requires a number 
of steps. The process starts with the development of a background environmental model. 
This model determines the path of acoustic propagation for the source to a receiver. The 
environmental model can be initialized with either historical data or in situ data obtained 
through traditional oceanography techniques. The acoustic data consists of a phase 
encoded sequence and timing information. Space-time signal processing is conducted to 
demodulate and decode the acoustic signal and provide spatial information for arrival 
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Figure 1.1: Acoustic tomography problem summary'. This thesis contribution is 
indicated by region in gray. 

identification. The environmental model influences the signal processing through 
beamforming and the expected arrival times. The array weights needed for a spatial 
matched filter for modal beamforming are provided by the environmental model. Arrival 
detection is performed in a temporal window around the expected arrival time provided 
by the environmental model. The difference between the actual arrival times as 
detemiined by the signal processing and arrival detection, and the expected arrival times 
is the information used by the inverse algorithm to determine the sound speed map. This 
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sound speed map may be used to improve the environmental model and an iterative 
solution may be determined. 

The objectives of this thesis are to; 

1. Develop space-time signal processing algorithms. 

2. Develop arrival detection, identification, and tracking algorithms. 

3. Assess beamformer performance. 

4. Process 13 hours of acoustic data from the 1992 Barents Sea Polar Front 
Experiment. 

The contribution of this thesis to the tom.ography problem is shown as the shaded 
region in Figure 1.1. The signal processing included plane wave beamforming, modal 
beamforming, and decoding of the signal. After the signal processing was complete, peak 
detection and travel time estimation were performed. All signal processing and peak 
detection algorithms were written in MATLAB™. MATLAB’^*’^ is a software package 
from the MathWorks Company. 

A description of the 1992 Barents Sea Polar Front Experiment is contained in 
Chapter II. The description includes the physical setup of the experiment and its goals. 
A general description of the oceanography is included for completeness. 

Chapter III summarizes the basic acoustic tomography theory. The discussion 
begins with propagation theory including both ray and mode theory. The discussion then 
shifts to the acoustic signal. Reasons for transmitting a maximal length phase encoded 
signal is discussed along with a description of these signals. Beamformering theory is 
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covered including both plane wave and modal. Finally all these pieces are tied together 
with a discussion of the tomographic inverse problem. 

Predictions for beamformer performance are covered in Chapter IV. These 
predictions included array beampattems and expected output. Synthetic data is used for 
later comparisons with experimental data. 

Chapter V discusses results from the processing of data. These results are compared 
with the predictions of Chapter IV. Results from inversions of data are included. 

Chapter VI summarizes the work completed and discusses lessons learned. Future 
work recommendations are given. 
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II. BARENTS SEA EXPERIMENT 



A. EXPERIMENT OBJECTIVES 

The Barents Sea Polar Front Experiment (BSPFEX) was conducted in August of 
1992. The experiment was a joint effort between the Naval Postgraduate School (NPS), 
Woods Hole Oceanographic Institution (WHOI), and the Science Applications 
International Corporation (SAIC). The principal investigators for the experiment are 
Professors Ching-Sang Chiu, James H. Miller, and Robert Bourke from NPS, Dr. James 
F. Lynch from WHOI, and Dr. Robin Muench from SAIC. The principal engineers for 
the development and deployment vertical hydrophone array system were Mr. Keith Von 
der Heydt and Mr. John Kemp from WHOI. The objectives of the experiment as outlined 
by the Barents Sea Polar Front Group [Ref. 5] were: 

1. Provide a detailed physical description of the polar front. 

2. Enhance the understanding of dynamics of the front, including frontogenesis and 
its influence on regional oceanographic processes. 

3. Assess the ability of acoustic tomography to define frontal and associated 
mesoscale features. 

4. Provide improved documentation of shallow water acoustic propagation in this 
region and the effect of the environment on acoustic Anti-Submarine Warfare 
(ASW) operations. 
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B. EXPERIMENT LOCATION 



The experiment was performed in the Barents Sea within a 80 x 90 km rectangle 
centered 100 km east of Bear Island as shown in Figure 2.1. The Barents Sea is bordered 
to the north by Spitzbergen and Franz Joseph Land, to the south by Scandinavia and 
Russia, to the east by Novaya Zemlya, and is open to the west. The experiment used 
both acoustic tomography and traditional oceanographic techniques to study the Barents 
Sea Polar Front [Ref. 6]. 




Figure 2.1: Experiment location. 
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The Barents Sea Polar Front is a stable feature which is caused by currents in the 
region [Ref 7]. These currents transport relatively cold fresh Arctic water and warm salty 
Atlantic waters to converge at the Barents Sea Polar Front. The characteristics of the 
different water masses are given in Table I and are shown in Figure 2.2. Figure 2.2 
depicts Lx)eng’s [Ref 8] model of the region. The abbreviations are: Arctic Water (AW), 
North Atlantic Water (NAW), Svalbard Bank Water (SBW), Barents Sea Water (BSW), 
and Bottom Water (BW). The experiment location is indicated by a rectangle in Figure 
2 . 2 . 



Table I: Barents Sea Water Masses [Ref. 6] 



Major Types 


T (”C) 


Salinity (psu) \ 


Arctic Water (AW) 


< .5 


34.2 - 34.0 1 


North Atlantic Water 
(NAW) 


-.5 to 2.0 


34.8 - 35.0 



i 

Local Variants 


T ("O 


Salinity (psu) 


Svalbard Bank Water 
(SBW) 


1.0 to 3.0 


< 34.4 


Bottom Water (BW) 


< -1.5 


> 35.0 


Barents Sea Water 
(BSW) 


-1.5 to 2.0 


34.7 - 35.0 


Polar Front Water (PW) 


-0.5 to 2.0 


34.8 - 35.0 
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Figure 2.2: Barents Sea water masses I Ref. 6|. Water masses as described In the 
text. 



The physical oceatiography of the region [days a large roie in controffing sound 
propagation |Rcf. 7). Ihc regioti is shaliow having an average depth of about 200 m. 
The bottom is rough in some areas and consists of a thin sediment overlying hard rock. 
The bathymetry and sound speed Is highly range dependent. The shallow water causes 
many surface and bottom interactions in the acoustic propagation. 
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The sound speed profile varies spatially [Ref. 6]. The sound speed profile measured 
at the vertical line array is shown in Figure 2.3. This sound speed profile is a typical 
Barents Sea profile south of the front for the summer. The sound speed profile gradient 
is basically negative from surface to bottom with a shallow mixed layer region. The 
negative gradient tends to bend rays toward the bottom. 



Sound Speed Profile near the vertical line array 




i^igure 2.3: Sound speed profile near the vertical line array. 
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C. EQUIPMENT DESCRIPTION 



The locations for moorings used in the experiment is shown in Figure 2.3 and the 
mooring locations given in Table II [Ref 9]. The moorings consisted of two 400 Hz 
transceivers (C and D), a 224 Hz source (A), and a vertical array of hydrophones (B). 
The characteristics of the sources is given in Table III. This thesis focus is on the signal 
processing of transmission from the 224 Hz source tothe vertical array. The bathymetry 
between the 224 Hz source and vertical line array are shown in Figure 2.4. 



Table II; Mooring Locations 



Mooring 


Latitude 


Longitude 


Depth (m) 


(B) SE(a) VLA 


74“ 19.1512’N 


23“ 33.1438’E 


275.0 


(B) SE(b) VLA 


74“ 19.1996’N 


23“ 32.2960’E 


275.0 


(A) NE 224 Hz 


74“ 37.5535’N 


23“ 24.3755’E 


142.0 


(D) NW 400 Hz 


74“ 32.9152’N 


21° 44.1043’E 


176.0 


(C) SW 400 Hz 


74“ 04.7337’N 


22“ 00.4605 ’E 


380.0 
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Table III: BSPFEX Source Characteristics 



Mooring 


(A) 


(C) 


(D) 


Frequency (Hz) 


224 


400 


400 


Source Level (dB) ref pPa 


183 


183 


183 


Bandwidth (Hz) 


16 


100 


100 


Duration (s) 


118.25 


132.86 


132.86 


Cycles/digit 


14 


4 


4 


Number of digits 


63 


511 


511 


Sequence length (s) 


30 


24(+2) 


24(+2) 


Repetition cycle 


2.5,7.5,... 


0,10,20,... 


5,15,25,.. 



The 224 Hz source mooring is depicted in Figure 2.5 [Ref 10]. The source 
consisted of quarter wavelength resonant tubes driven by a piezoelectric transducer. A 
phase encoded signal (which is discussed in detail in Chapter III), was transmitted from 
the 224 Hz source and received on the vertical line array. 
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Figure 2.4: Mooring locations for tomography experiment. 
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3D Representation of Region of Study 




23.7 



74.35 



Latitude 



23.1 



23.2 



Longitude 



figure 2.5: Three dimensional representation of transmission path for the tomography 
experiment. 
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The vertical line array is shown in Figure 2.7 [Ref. 11]. The array consisted of 16 
hydrophones spaced 10 meters apart. The hydrophones have a sensitivity of -160 dB re 
V/pPa. The mooring was designed to acoustically decouple the surface buoy from the 
array and minimize the effect of surface waves on mooring motion. The surface buoy 
contained an RF Ethernet system [Ref. 11] which transmitted the digitized data to the 
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ship. On the ship the data was stored on tape for later analysis. Some hydrophones 
experienced a loss in sensitivity due to an electrical problem. These hydrophones were 
later identified and were compensated for in the analysis of the data. 




"igure 2.7: Diagram of the WHOI-NPS vertical line array used for the Barents Sea 
Polar Front Experiment [Ref. II]. 
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III. ACOUSTIC TOMOGRAPHY THEORY 



A. TOMOGRAPHY PROBLEM OVERVIEW 

The acoustic tomography problem consists of two parts, the forward problem and 
the inverse problem [Ref 2]. The forward problem establishes the physical relationship 
between data and the unknown structure. The inverse problem consists of estimaring the 
unknown structure using data collected. 

For acoustic tomography, the forward problem consists of a detailed analysis of the 
expected acoustic propagation for the area to be studied. This modeling of the 
environment determines the path of propagation and the expected arrival times for 
acoustic signals. Figure 3. 1 shows a simplified view of the forward problem. Note that 
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figure 3.1: Simplified picture of shallow water tomography problem. There are 
many paths between source and receivers. 
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two of the paths have the same path length. Even with a background field that has a 



range dependent sound speed profile, two arrivals may occur too close together in time 
to distinquish. For the worst case, an isovelocity sound speed profile, there would be two 
arrivals at the same time. To distinguish these different arrivals, upward and downward 
propagating sound waves must be distinguished. Therefore in shallow water tomography, 
an array of hydrophones is often necessary for to separate the arrivals at different angles 
[Refs.3, 4]. 

The sound velocity field can be thought of as the sum of a background sound field 
with a sound speed perturbation [Ref 12], 



If 5c is small compared to c„, then using a binomial expansion and substituting Equation 
3.1 in to Equation 3.2 yields. 



c{x,y,z,t) =c^(x,y,z) + 5c(x,y, z, t) . 



(3.1) 



The travel time from the source to receiver along some path is given by. 




(3.2) 




(3.3) 



Equation 3.3 can be approximated as 
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(3.4) 






6c 



cl{x,y,z) 



ds. 



The arrival time can now be broken down into an expected arrival time Tjo and a 
perturbation 6Tj as. 



Xj = + 6x^. 



(3.5) 



The arrival time perturbation is then given by 



5x, = - f 

path! cl{x,y,z) 



(3.6) 



The inverse problem is then to determine he from the arrival time perturbation 6x. The 
determination of the time perturbation is then dependent on the solution of the forward 
problem which is given by. 



X 



i,o 



I 

path i 



1 

c^{x,y, z) 



ds. 



(3.7) 



The data for the inverse problem is 5Xj as indicated by Equation 3.6. The data is the 
difference between the measured arrival time of acoustic signals and the model predicted 
arrival times. The model arrival times are used in the detection and identification of 
actual arrivals. The model arrival times focus in where to look for the expected arrivals. 
Two types of models were used, a ray (plane wave-like) propagation model and modal 
propagation model. 
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B. FORWARD PROBLEM 



1. Plane Wave Propagation 

Ray theory is the result of an approximate solution to the linearized wave 
equation [Ref. 12]. The wave equation is given by 

V^P - " -47t6(f - (3.8) 



where is the source location and the sound speed and pressure can vary 

spatially and temporally. A solution to the wave equation is 



p{x,y, z, t) =A{x,y, z) e 



ju(t-ru,y, 2)/Ci) 

t 



(3.9) 



where A is the pressure amplitude, c, is a constant phase speed, and T has units of length. 
The direction of travel for the wave is given by VT. Substituting Equation 3.9 in to 3.8 
yields: 



^-(i!>)^vr*vr + ( 2 — •VT+v^D = o. (3.io) 

Ac C A 



If the following approximations are made: 



sfa 



fe)’ ^ 



fV2ri<-“ 

C 






(3.11) 



then the wave equation can be approximated by: 
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W»W = , n{x,y,z) = - .- -2 . ? . — (3.12) 

c{x,y, z) 

where n(x,y,z) if the index of refraction and c„ is an arbitrary reference sound speed. This 
is known as the Eikonal equation. The approximations physically mean: 

1. The amplitude must change slowly compared to a wavelength. 

2. The speed of sound must change slowly compared to a wavelength. 

The behavior of VT can be found from 

^ (VT) = Vrt. (3.13) 

ds 

The model used for ray tracing was the NO A A Hamiltonian Acoustic 
Raytracing Program for the Ocean (HARPO) [Ref. 14]. Work on the ray forward 
problem was done by LT John Mykyta, USN [Ref. 6]. Programs external to HARPO 
were used to perform eigenray searches and account for losses from spreading or 
boundary interactions. Spreading losses, surface scattering losses, and bottom losses were 
incorporated after the raytracing was completed [Ref 6]. 

2. Modal Propagation 

Equation 3.8 can also be solved for normal modes. Expressing Equation 3.8 
in cylindrical coordinates yields [Ref. 15] 

- -|6(r)6(z-z,)e-^ (3.14) 

With separation of variables, the solution to Equation 3.14 can be expressed as 
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(3.15) 



p(r, z) = Y, 



n=l 



R„{r) Z^(z, r) 
yfr 



where Z, represents normal modes and is a solution to a Helmholtz equation [Ref. 15]: 

Z„ = Q, 

where !^(r,z) is the wave number and is the eigenvalue. also represents an 
orthonormal set which satisfies 





-F 


^ {r , z) _ 


|p^(z) dz) 




p^iz) pjz) J 



/ 



^I^dz = 6. 



Po 



(3.17) 



The radial component satisfies [Ref. 15]: 



dr^ 






18) 



Neglecting ll(4r^) term and assuming adiabatic the solution becomes: 



z 

(3.19) 

yf^ 

For normal mode modeling, mode coupling must be considered for significantly range 
dependent problems [Ref. 3]. 



p(r, z) = 



C. SPACE-TIME SIGNAL PROCESSING 
1. Temporal Signal Processing 
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In acoustic tomography, arrival times must be measured to within a few micro 
seconds. In order to achieve this accuracies, impulse-like signals are desirable with high 
signal-to-noise ratios [Refs. 4,12]. Impulses are ideal for mathematieal modeling but are 
physically difficult to achieve with finite bandwidth sources whose peak power is limited. 
One solution is to use pulse compression techniques. The idea is to spread the energy out 
in the time domain efficiently using the available bandwidth of the source. At the 
reeeiver, the energy is combined or eompressed to an impulse-like arrival. Coding 
schemes such as maximal-length sequences (m-sequences) have been used with success 
in a number of recent tomography experiments [Refs. 4,12]. 

M-sequences can be constructed from primitive polynomials [Ref. 16]. The 
signal for the 1992 Barents Sea tomography experiment [Ref. 4] was based on a sixth 
degree primitive polynomial given by 

g{D) = + £) + 1 . (3.20) 

The output of the shift register (Figure 3.1) taken at any of the delays is a unique, 
periodic, binary sequence of length N^, = 2“ - 1 = 63 before repeating. The actual 
sequence used for the Barents Sea tomography experiment is shown in Figure 3.2. 

M-sequences have a number of interesting properties which are outlined in a 
text book by Ziemer and Peterson [Ref 16]. An important property for the purposes of 
tomography is the autocorrelation of the sequence. If the binary sequence (0,1) is mapped 
to a (1,-1) sequence, then the normalized autocorrelation is shown in Figure 3.3. Note 
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that the autocorrelation consists of a triangle with a width of two digits and height of one 
and is everywhere else HN„. 

Birdsall and Metzger (referred to as BM) [Ref. 17] have pioneered the use of 
m-sequences in acoustic tomography. To build a BM signal, the phase is encoded with 
the m-sequence. The phase-encoded signal is 

s{t) = Acos{2%f^t m{t)d^) (3.21) 

where ^ is the carrier frequency, /I is a constant amplitude, m(t) is the m-sequence to be 
encoded, and 0^ is a constant. The constant 0„ can be chosen to maximize the relative 
power between the carrier and sidebands. The 0^ which maximizes the output signal-to- 
noise ratio is given by [Ref 17] 

0^ = tan-i(^ (3.22) 

which for the Barents Sea experiment was 82.8192 degrees [Ref. 4]. The frequency 
spectrum for the signal is shown in Figure 3.4. The complex envelope of the BM signal 
is [Ref. 4] 

s{t) = (3-23) 

The important feature of the autocorrelation of the complex envelope of the BM signal 
is the lack of sidelobes, making arrival time estimation of such a signal as accurate as 
possible. This is seen by computing the periodic autocoirelation as 
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(t) = f 






dx , 



(3.24) 



T=N„d 



where the integration is over one period (T = N„d) of the m-sequence. BM show that this 
autocorrelation is 0 outside the main lobe as shown in Figure 3.5, i.e. 



p (r\ = if X = 0 



(3.25) 



Equation 3.24 can be converted to the frequency domain as 



R^,{f) = S{f)S^{f) 



(3.26) 



where S(f) is the fft of s(t). Since beamforming was conducted in the frequency domain, 
temporal filtering was easily incorporated with Equation 3.26. Prior to processing, a 
replica of the signal s(t) was converted through an fft to S(f). The complex conjugate of 
S(f) was then used for circular convolution in the frequency domain. 
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Figure 3.1: Delays for generating m-sequence (modulo 2 sum). 
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M-sequence for 224 Hz source 




figure 3.2: M-sequence code for 224 Hz source. 
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Figure 3.3: Autocorrelation property of m-sequence. 
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Frequency Spectrum for m-sequence 




figure 3.4: Coded signal frequency spectrum. 



Temporal resolution kernel of BM processing for 62,5 msec digit 




Figure 3.5: Temporal resolution kernel of BM processing for 62.5 msec digit [Ref. 4]. 
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2. Spatial Signal Processing 



a. Plane Wave Beamforming 



The resolution and identification of individual rays can often be 



accomplished by plane-wave beamforming the outputs of individual hydrophones on a 
vertical array [Ref. 4]. The outputs of the hydrophones can be phase shifted to maximize 
the sensitivity in a particular direction. The complex envelope of the n* ray, which is 
assumed to be a local plane wave acoustic signal with m-sequence phase modulation, is 
given by [Ref. 4] 



where i! = (u,v,w) are the direction cosines with respect to the (x,y,z) axes and are given 
by 



in terms of the spherical angles 0 and \)/. The Fourier transform of the complex envelope 
of the n'*' ray’s pressure field is 



where Sn,(f) is the Fourier transform of the complex envelope of the m-sequence phase 
encoded signal given by 




u = sin6cost|; 
V = sinGsintj; 
W - COS0 



(3.28) 




(3.29) 
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S^(f) = . (3.30) 

The venical array beamforming and associated ra-sequence processing can be thought of 



as a spatial-temporal matched filter which as a transfer function is given by 

B^^y{f,x) = ^ 6 ( 2 - 2 ^) 6 {x-x^) 5 (y-y^) S*(f) (3-31) 

where H is the number of hydrophone, z^, are the hydrophone depths, x, and y, are the 
coordinates of the vertical array, and Tj, is an adjustable time delay applied to each 
hydrophone. The complex envelope of the output of any general beamformer is given by 



Y{f) ^ Jp{f,x)‘ { f, x) dV^ (3.32) 

where is the volume of the aperture or array. For the Barents Sea Tomography 
Experiment we have 



Y{f) 
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where a local coordinate system is assumed for the vertical array. If the time delays are 
adjusted for 



t. = -iicose,, (3.34) 

^ c 

we have a spatial-temporal matched plane-wave filter for the n* ray. If the acoustic field 
is well modeled, different delays may be included such as wavefront curvature. 
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The uncertainty in the estimate for travel time depends on the signal to 
noise ratio [Ref. 4]. The standard deviation of the travel time estimation error is well 
modeled by 



= 



d 

yfENR' 



(3.35) 



where SNR is the signal-to-noise ratio and d is the digit width of the maximal length 
sequence modulating the phase of the transmitted signal. Equation 3.35 is valid for SNR 
greater than 10. 

b. Mode Beamforming 

The resolution and identification of individual normal modes can be 
accomplished by beamforming the outputs of individual hydrophones on a vertical array 
[Ref. 3,4]. Instead of time delaying the outputs by an appropriate amount (or phase 
shifting in the frequency domain) to plane- wave beamform for a ray, the amplitude of the 
mode at each of the hydrophones across the band of frequencies the signal occupies is 
matched. The temporal Fourier transform of the complex envelope of the n* mode’s 
pressure field is 



Pjf,x) = (3.36) 

where Zffpc,y,z) is the nth local mode function which satisfies orthogonality in depth at 
all (x,y) and SJf) is the Fourier transform of the complex envelope of the m-sequence 
phase-encoded signal as in the ray case. The vertical array and m-sequence processing 
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can be thought of as the spatial-temporal matched filter and has a transfer function given 



by 



H 



= g6(z-z^)6(x)6(y)Z^(f,x,y, (3-37) 

where H is the number of hydrophones and are the hydrophone depths. The complex 
envelope of the output of the mode beamformer is 



Y{f) = jp{f,x) { f, X) dv^ 



(3.38) 



For the Barents Sea Tomography Experiment the modal beamformer can be written as 
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Y{f) = J]a„Zjif,zjZjz)5^(if)5;(if)e 



-j2nft^ 



(3.39) 



h=l 



The assumption that the modes are orthogonal over the array aperture 
leads to errors since the array only covers a portion of the water column. Shang [Ref. 
15] suggests a least squares approach to modal beamforming. The least square problem 
is to find weights A(f) such that the estimated field matches the measured field. The 
derivation for these weights A(f) are shown here for completeness [Ref 18]. Assuming 
that the estimator A(f) is linearly related to the observations yields 

A = Bd, (3-40) 

where d is the data, and 5 is a matrix to be determined. The error is given by; 

e=i-A. (3.41) 



The covariance of the error is 
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Cg = E[ee'^] . 



(3.42) 



The constant B is to be chosen to minimize Q. The covariance of the error 
written as 



= E[ {Bd-A) {Bd-A) n 
= E[ {Bd-A) id^B^-A^) ] 

= E[Bdd^B^ - Ad^B^ - BdA^ + AA ^ 
= BC^^ - - BCla + 



or equivalently as 



Ce= iB - C^{B - ^ + c^. 



The minimization of Q gives , 



B = 






Substituting Equation 3.45 in to Equation 3.40 yields 



A = 



where the covariance A and d, and Cj is the covariance of the data d. 
Assuming that the data can be represented as 



can be 



(3.43) 



(3.44) 



(3.45) 



(3.46) 
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d = M{z,f)A{f) +v„{z,f) 



( 3 . 47 ) 



where M(zf) is a matrix of mode weights at a hydrophones, A(f) is a matrix of 
amplitudes, and vjzf) is noise associated with the observation at a hydrophone. The 
modes weights matrix is written as 



M(z,f) = 



'm. 


/O 


^,0 


^ 3,0 • 




M, 


A 


^ 2,1 


^ 3,1 • 




Ml, 


15 


^2,15 


^ 3,15 • 


• • 15 , 



( 3 . 48 ) 



where the first subscript indicates mode number and the second indicates hydrophone 
number. The amplitude matrix is written as 



■^1 

^2 



A{f) = ' 



( 3 . 49 ) 



where the subscript indicates mode number. The least-squares problem is then to find the 
amplitude A(f) where 
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A = A(d) . 



(3.50) 



The solution to 3.50 is given by equation 3.42. The remaining portion of the problem is 
to find and Cj. C^d can be found as follows: 



= E[A{MA+v)'^ 

= E[A{A'^M^+v'^) ] 
= E[AA'^M'^ + Av'^] 



Cd is found as foUows: 

= E[{MA + v) {MA + v) n 
= E[{MA v) {A + v^) ] 

= E[MAA'^M'^ + vA'^M'^ + MAv"^ + w"^] 
= MC^M^ + CavM^ + + C^. 



If the noise v and modal amplitudes are assumed to be uncorrelated, then Equation 3.52 
simplifies to 



= MCj^M^ + C^. (3.53) 

The solution can now be written by substituting Equations 3.53 and 3.51 into Equation 
3.42 which yields 



A = {C^M + + Cj-"d. (3.54) 

Equation 3.54 can be manipulated into a more familiar form of 

A = (CX^ + '^Af^c;^d, (3.55) 

where Q represents any a priori information, Cy is the covariance of the noise, and C, 
represents the covariance of the error in the estimate. 
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The acoustic field at the array is described in a least-squares sense as 



(9 =Y^M{z, f)A{f) , 

m 



(3.56) 



where A(f) was estimated as shown in Equation 3.55. The best beamformer in a least- 
squares sense is then described by 

Y{f) = a„Zjf, z^) Zjz) Sjf) s:{f)A{f) . (3-57) 

/i=i 



As with Equation 3.39, the beamformer is summed over the 16 hydrophones. 
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IV. PERFORMANCE PREDICTIONS 



A. PLANE WAVE PREDICTIONS 

The hydrophone array used for the Barents Sea tomography experiment consisted 
of 16 elements separated by 10 meters as shown in Figure 2.7. The beampattem for the 
array at 224 Hz is shown in Figure 4.1. Note the grating lobes are at plus and minus 40 
degrees. These are caused by the spatial undersampling by the array. The 224 Hz 
frequency has a wavelength of approximately 6.5 meters and would require a hydrophone 
spacing less than 3.25 meters to prevent grating lobes. The grating lobes are acceptable 
because the arrivals are separated in both time and space. Thus energy coming in through 
a grating lobe will not interfere temporally with the beamformed arrival. A rectangular 
window was chosen for the array to achieve the narrowest beamwidth possible for arrival 
identification. The -3 dB beamwidth for the array is 2 degrees at 224 Hz. 

Chapter II showed the physical setup for the tomography experiment. The sound 
speed profile (Figure 2.3) is basically a negative gradient. The effect of this environment 
for plane waves was modeled by John Mykyta using HARPO and other codes [Ref 6,14]. 
A plane wave assumption leads to errors since in a non-isovelocity environment the 
arrivals will not be plane waves. To show the effect of plane wave beamforming when 
the incoming wave is not a plane wave, a one millisecond delay was allowed over the 
array. The result is shown in Figure 4.2. The effect is that sensitivity is reduced by one 
dB and the sidelobes are higher. If the environment is known well enough, these 
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wavefront curvatures can be taken in to account as pointed out in Equation 3.19. Several 
approaches are available. The corrections for wavefront curvature can be done using only 
the background environmental model or after each inversion, produce a new 
environmental model for wavefront correction. 




Figure 4.1: Plane wave beampattem for Barents Sea 



VLA at 224 Hz. 



The actual amount of wavefront curvature varies with arrival angle. An example 
of the predicted wavefront curvature is shown in Figure 4.3. The wavefront curvature 
predictions were done by John Mykyta using HARPO [Ref. 6]. These predictions were 
made while determining the arrival structure at the array. 

The difference between the predicted arrivals and the actual arrivals is used to find 
the sound speed perturbations. The predicted arrival structure is shown in Figure 4.4. 
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While the arrival structure appears to be smooth, it does not take in to account the 
beampattem of the array. A more representative picture of the arrival structure is shown 
in Figure 4.5. Here the arrivals within a one degree bin are shown. For low angles, there 
is little separation between the arrivals, but for the higher angles (greater than plus or 
minus eight degrees), there is sufficient temporal separation to resolve more than one 
arrival. Figure 4.6 shows the predicted beamformed arrivals. Thus for higher angles 
multiple arrivals should be seen. Arrivals are indicated by the higher (red) amplitude 
regions. This figure was produced by beamforming synthetic data. Multiple arrivals are 
evident in higher angles. 



Normalized Beampattem using Boxcar Window 




Figure 4.2: Wavefront curvature degradation prediction. 
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Example of HARPO model prediction 




Figure 4.3: Example of predicted wavefront curvature. 
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Predicted plane wave arrival structure 




Figure 4.4: Predicted plane wave arrival structure [Ref. 6]. 
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Predicted arrival structure in one degree bins 




Figure 4.5: Predicted plane arrival structure in one degree bins. 
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Predicted Beamformed Ray Arrival Structure 




Travel time (sec) 



Figure 4.6: Predicted plane wave beamformed arrival structure [Ref. 3]. 
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B. MODAL BEAMFORMER PREDICTIONS 



Modal beamformer performance prediction is more difficult. The modal 
beamformer depends upon using a spatial filter to discriminate modes. The weights for 
the spatial filter depend upon knowing the mode shapes. Mode shapes at the array are 
dependent upon the sound speed profile and the boundary conditions. Performance will 
depend upon how well the environment is known at the array and how much the 
environment changes. 

The modal beamformer is described by Equation 3.21. The performance of the 
beamformer depends upon the modes being orthogonal as described in Equation 3.10. 
A measure of the cross talk at a given frequency can be determined by approximating 
Equation 3.10 over the extent of the array aperture as: 



r 



n,m 



Po ic = 0 



Znik)Zjk) 



(4.1) 



where ^ is the hydrophone depth, n and m are mode numbers, and the result is interpreted 
as the fraction of energy from mode n when beamforming for mode m. Ideally the result 
would be zero unless n=m in which case the result would be one. Using Equation 4.1 
at the center frequency of 224 Hz, Figure 4.7 was generated. This figure shows that 
modes one through three are well resolved. Modes one through three all having turning 
depths within the array so are expected to perform well [Ref. 9]. For higher modes, 
crosstalk is expected. 

As with the plane wave beamformer, temporal separation of the modes also will 
assists in identification. Figure 4.8 shows the expected modal arrival structure. Figure 
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4.9 shows the ideal modal beamformed expected arrival structure. Figure 4.9 would be 
achieved if the vertical line array covered the entire water column. Figtire 4.10 shows 
the effect of the beamformer. Comparing Figure 4.9 and Figure 4.10 we see the arrivals 
are visible and the timing remains intact but the cross talk is evident with more than one 
arrival visible. 




Figure 4.7: Predicted modal beamformer crosstalk. 
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Predicted modal arrival structure 




Figure 4.8: Predicted modal arrival structure [Ref 3]. 
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Figure 4.9: Ideal modal arrival structure. The "perfect" beamformer would have no 
cross talk [Ref 3]. 
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Figure 4.10: Beamformed synthetic data. This figure shows the expected cross talk 
between modes. Temporal separation makes arrival identification possible [Ref. 3]. 
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V. RESULTS 



A. PLANE WAVE 

Plane wave beamforming was performed on 13 hours 20 minutes of acoustic data. 
For each five minutes of data, there were 30 transmissions of 3.9375 seconds. Of these 
30 transmissions, the middle 28 were analyzed. The first and last transmissions were not 
used because the pulse compression of the signal was done in the frequency domain thus 
requiring a periodic signal. Using circular convolution, the first and last sequence may 
not have acoustic data over the time analyzed. After beamforming each of the 28 seconds 
of data from -19 to +18 degrees, the 28 sequences were coherently averaged prior to peak 
detection, identification, and tracking. 

Figure 5.1 shows a sample of decoded signals on hydrophone seven for one set of 
30 transmissions. Only the middle 28 arrivals are shown. The arrivals show a large peak 
in the beginning which is from low angle arrivals. Beamforming is necessary to 
distinguish different individual arrivals in the initial large peak and to bring out later 
arrivals. After beamforming, the 28 arrivals are coherently averaged. This figure shows 
the consistency in the signal over the different transmissions. 

Figure 5.2 shows an example of plane wave beamformed output. The figure is 
windowed for a region of time which the arrivals are expected. The actual arrivals are 
indicated by the red regions. Note that multiple arrivals are apparent as was predicted 
in the forward model. It is the difference between the arrivals shown in Figure 5.2 and 
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Figure 4.6 which provide the information for the inverse algorithm. Figure 5.3 shows an 
example of actual detected peaks and modeled peaks. Using information as shown in 
Figure 4.3, sound speed maps were produced for each five minutes [Ref. 3]. An example 
of a sound speed map is shown in Figure 5.4. 



Decoded signals for 28 transmissions on hydrophone 7 




figure 5.1: Example of decoded acoustic signal. 



50 



229:04:37 UTC, Beamformed Ray Arrivals 




Figure 5.2: Example of plane wave beamformed output. 
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Predicted and actual arrivals for day 0229 04:37 




Figure 5,3: Example of model predicted arrival times with detected arrival times for 
plane wave beamformer. Predicted ’o’, actual ’x’. Predictions as in Ref 15. 



Figure 5.4 clearly shows the Barents Sea Polar Front. The Front has a fine vertical 
structure. The warm and salty Atlantic water is to the south and is indicated on the map 
by the red and yellow regions. From the processed acoustic data, 160 frames were 
produced similar to Figure 5.4 [Ref. 3]. These frames were then transferred to video tape 
producing a movie of the frontal dynamics over the tidal period. A summary of the data 
used to produce this video is show in Figure 5.5. 
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figure 5.4: Example of sound speed map generated by the inverse algorithm 
developed by Chiu [Ref. 3]. 
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beamforming. This data was used to produce a sound speed movie. 
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Wavefront correction was attempted for the plane wave beamforming but no 
significant improvement was noted. The wavefront correction was determined using 
model wavefront curvature for the background sound speed field [Ref 6]. The model 
prediction for the wavefront curvature are difficult because of some limitations of the 
model. Small local variations in the bathymetry results in arrivals which are not smooth. 
This is because HARPO treats the rays as having no radius [Ref 14]. Thus the rays do 
not ensonify a region of the bottom but instead reflect off a point and are not influenced 
by the region. 

B. MODAL BEAMFORMING RESULTS 

Modal beamforming was performed on the 13 hours and 20 minutes of data. Two 
types of beamforming was done. The first modal beamforming was an outer product 
performed as indicated in Equation 3.29 with its associated beampattem shown in Figure 
4.7. This beamforming yielded identifiable arrivals. Figure 5.6 shows an example of 
modal beamformed output. Note the similarities between Figures 5.6 and 4.9. The 
arrivals can be identified as the red regions near the center of the windows. Windowed 
regions are shown to make arrivals easier to identify. Evidence of the crosstalk between 
modes can be seen in both Figures 4.9 and 5.6. Figure 5.7 shows an example model 
predicted arrival times and beamformed arrival times. The difference between these is 
the information that can be used for the inverse algorithm to determine spatial sound field. 
A summary of travel modal travel time results is shown in Figure 5.8. 
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A least-squares approach was also used for modal beamforming. Using Equation 
3.46 weights were determined using least-squares for the first 14 modes. Higher modes 
were not analyzed because the solution became ill conditioned for more than 14 modes. 
While this method produced results which were consistent, the outer product beamformer 
produced results which could be used to identify and track arrivals for 30 modes. 



229:04:37 UTC, Beamformed Mode Arrival 




Travel time (sec) 



Figure 5.6: Example of modal beamformed output. 
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Example of predicted and actual modal arrivals 
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Figure 5.7: Example of modal travel time prediction and actual detected travel times. 
Model predictions indicated by ’o’, actual travel time indicated by ’x’. Predictions 
from [Ref. 3]. 



Modal Travel Time Data 




Figure 5.8: Summary of modal arrival travel times. This data can inverted to 
produce sound speed maps. 
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VI. SUMMARY AND FUTURE WORK 



A. SIGNAL PROCESSING SUMMARY 

Temporal and spatial signal processing algorithms were developed and applied to 
acoustic data acquired in the Barents Sea Polar Front Experiment in 1992. Over 13 hours 
of acoustic data was processed for plane waves without wavefront corrections and over 
eight hours was processed for normal modes. An arrival detection, identification, and 
tracking algorithm was developed and applied to both modal and plane wave beamformed 
data. Results without wavefront correction were predicted to be adequate to resolve 
arrivals. The output of the plane wave beamformer without wavefront correction yielded 
adequate results to develop tomographic images. 

Wavefront correction was attempted but while the ray model produced accurate 
results for arrival estimation, model limitations prevented accurate wavefront information. 
Tracing rays at finer launch angles may yield sufficient data to perform wavefront 
correction. With more data smoothing could be done to produce reasonable wavefront 
corrections. There are indications of a great deal of variability in the sound speed map. 
The background model may not provide accurate enough predictions over a long period 
of time to correct for wavefront curvature. 

Modal beamforming was conducted using a dot-product and least-squares approach. 
The dot-product beamformer produced adequate results to determine arrival information 
for 30 modes. These results should be adequate for future inversion. The least-squares 
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beamformer was limited to the first 14 modes because the solution became ill conditioned 



for more. An array with more elements should allow a least-squares solution for more 
modes. The least-squares approach should produce a more statistically reliable 
beamformer but more array elements are needed for a better solution. 

B. FUTURE WORK 

The algorithms are complete for plane wave and modal beamforming. The data 
which was processed was chosen because it was the largest continuous block of time. 
There are other blocks of time which could be processed. After completion of processing 
a linear inteipolation could be done to extend the period of time covered by the movie. 

Wavefront correction should be analyzed further. The first cut would be to perform 
ray tracing using HARPO with smaller differences in launch angle, beginning with 0.001 
degrees between launch angles. Then analyze the data for improvement. A second cut 
would be to perform ray tracing using the sound speed maps which have been previously 
generated. These sound speed fields should produce more accurate ray tracing results. 

Perform a detailed study on least-squares beamformer performance. Begin with 
synthetic data with a known noise field added. Then determine if the least-squares 
beamformer produces better results as measured by travel time error estimation. 

The beamformer code is currently written to be able to run on computers with 
relativly small amounts of RAM. With access to better computers, the beamformer code 
can be optimized to run much faster. Program speed was limited by read/writes to hard 
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disks for temporary storage. Eliminating unnecessary disk writes should speed the 
execution time of the program significantly. 

C. RECOMMENDATIONS 

The greatest impact on modal beamforming performance is the sound channel 
coverage. A longer array with more elements would improve the modal beamformer. 
The modal beamformer also depends upon accurately knowing the mode shape at the 
hydrophones. This requires knowing where the hydrophones are in depth. Tilt 
information was available for the array and a linear tilt model was used. The vertical 
array likely has curvature from currents and tilt. Better array shape information would 
improve processing through accurate knowledge of hydrophone position. 
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